%% Aim
% Preprocess (left) pupil data and generate a single pupil.csv

%% setup
tic;
clc;
clear;
root_path = pwd; 
addpath(genpath(root_path));

nrole=2;
nblock=2;
ntrial=100;
nprice=10;
nlottery=10;
  
%% load data
sub_data=readtable([root_path,'\data\subject\subject.csv']);
eye_data=readtable([root_path,'\data\eye\eye.csv']);
nsub=height(sub_data);

%% Clean data
eye_data.PupilLeft(eye_data.PupilLeft<=0)=NaN;

%% Generate 0) Remove extreme data, 1) interporlated, 2) filtered, and 3) zscored pupil for each subject 

for isub=1:nsub
    for irole=1:nrole
        
        edata=eye_data(eye_data.Subject==isub & eye_data.Role==irole,:);
    
        % Removing data of extreme sample-to-sample velocity
        % Kret, M. E., & Sjak-Shie, E. E. (2018). Preprocessing pupil size 
        % data: Guidelines and code. Behavior research methods, 1-7.
        edata.PupilDev=max(abs([NaN;diff(edata.PupilLeft)])./abs([NaN;diff(edata.TimeSignal)]),abs([diff(edata.PupilLeft);NaN])./abs([diff(edata.TimeSignal);NaN]));
        medPupilDev=nanmedian(edata.PupilDev);
        MAD=nanmedian(abs(edata.PupilDev-medPupilDev));
        TMAD=medPupilDev+2*MAD;
        edata.PupilValid=edata.PupilLeft;
        edata.PupilValid(edata.PupilDev>TMAD,:)=NaN;
 
        % Removing "island"
        edata.Isnan=isnan(edata.PupilValid);
        diffs = diff(edata.Isnan);
        idx=find(diffs~=0);
        idxnan=find(diffs==1); % change from num to NaN
        idxnum=find(diffs==-1); % change from NaN to num
        idxmax=height(edata);  
        
        for i=find(idx==idxnum(1)):2:find(idx==idxnan(end))
            pre_data=edata(edata.TimeSignal>edata.TimeSignal(idx(i)+1)-75 & edata.TimeSignal<edata.TimeSignal(idx(i)+1),:);
            post_data=edata(edata.TimeSignal<edata.TimeSignal(idx(i+1))+75 & edata.TimeSignal>edata.TimeSignal(idx(i+1)),:);
            pre_post_data=[pre_data;post_data];
            dur=edata.TimeSignal(idx(i+1))-edata.TimeSignal(idx(i)+1);            
            if dur<25 || (dur<50 && sum(~isnan(pre_post_data.PupilValid))==0)
               edata.PupilValid((idx(i)+1):idx(i+1))=NaN;
            end      
        end        

        % Resample to 2 ms interval by interpolation, default linear
        pupil_data=table;       
        [pupil_data.PupilValid, pupil_data.TimeSignal] = resample(edata.PupilValid,edata.TimeSignal,0.5,'linear');
        pupil_data.Subject(:)=isub;
        pupil_data.Block(:)=(2-(sub_data.Role1(isub)==irole));        
        pupil_data.Role(:)=irole;
        pupil_data.PricePosition(:)=sub_data.PricePosition(isub);
        
        % Relabel stimulus event
        pupil_data.StimulusName(:)=NaN;
        diffs=diff(edata.StimulusName);
        idx=find(diffs~=0);
        pupil_data.StimulusName(pupil_data.TimeSignal<=edata.TimeSignal(idx(1)),:)=edata.StimulusName(idx(1));
        for i=2:length(idx)
            pupil_data.StimulusName(pupil_data.TimeSignal<=edata.TimeSignal(idx(i)) & pupil_data.TimeSignal>edata.TimeSignal(idx(i-1)),:)=edata.StimulusName(idx(i));
        end
        pupil_data.StimulusName(pupil_data.TimeSignal>edata.TimeSignal(idx(end)),:)=edata.StimulusName(idx(end)+1);             
        
        % Filter with 3th order Butterworth, low pass, 4HZ
        [b,a] = butter(3, 4/(500/2));
        pupil_data.PupilFilter = filter(b,a,pupil_data.PupilValid);
        
        % Zscore
        pupil_data.Pupilz=zscore(pupil_data.PupilFilter);     
        
        % Save preprocessed data
        writetable(pupil_data, [root_path, '\data\pupil\s',num2str(isub),'_r',num2str(irole),'_pupil.csv']);
    
    end
end

%% generate the concatenated .csv file

pupil_data=[];
for isub = 1:nsub
    for irole = 1:nrole
        pupil_data = [pupil_data; readtable([root_path,'\data\pupil\s', int2str(isub), '_r', int2str(irole), '_pupil.csv'])];
    end
end
writetable(pupil_data, [root_path, '\data\pupil\pupil.csv']);

%%
toc;
